home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-02 / tsptp.zip / FBENCH.PAS < prev    next >
Pascal/Delphi Source File  |  1993-04-09  |  13KB  |  379 lines
  1. (******************************************************************************)
  2. (*                                FBENCH.PAS                                  *)
  3. (*                                                                            *)
  4. (*  This is a complete optical design raytracing algorithm, stripped of its   *)
  5. (*  user interface and translated into "standard" Pascal.  It not only        *)
  6. (*  determines execution speed on an extremely floating point (including trig *)
  7. (*  function) intensive real-world application, it checks accuracy on an      *)
  8. (*  algorithm that is exquisitely sensitive to errors.  The performance of    *)
  9. (*  this program is typically far more sensitive to changes in the efficiency *)
  10. (*  of the trigonometric library routines than the average floating point     *)
  11. (*  program.  This program has been demonstrated to be invariant under changes*)
  12. (*  in floating point format, as long as the format is a recognised double    *)
  13. (*  precision format.  If you encounter errors, please remember that they are *)
  14. (*  just as likely to be in the floating point editing library or the         *)
  15. (*  trigonometric libraries as in the low level operator code.                *)
  16. (*                                                                            *)
  17. (******************************************************************************)
  18.  
  19. PROGRAM FBench(Output);
  20.  
  21. (******************************************************************************)
  22. (*                                TIMING                                      *)
  23. (******************************************************************************)
  24.  
  25. (*$IFNDEF TopSpeed *)
  26.  (*%F TRUE   *** Compile for Turbo Pascal ***)
  27.   USES TPBench;
  28.  (*%E*)
  29. (*$ELSE     *** Compile for TopSpeed Pascal ***)
  30.   IMPORT TSBench *;
  31. (*$ENDIF *)
  32.  
  33. (******************************************************************************)
  34.  
  35.   CONST
  36.     MaxSurfaces = 10;
  37.  
  38.   VAR
  39.     CurrentSurfaces : BmInt;
  40.     Paraxial        : BOOLEAN;
  41.     ClearAperture,
  42.     AberLSpher,
  43.     AberOsc,
  44.     AberLChrom,
  45.     MaxLSpher,
  46.     MaxOsc,
  47.     MaxLChrom,
  48.     RadiusOfCurvature,
  49.     ObjDist,
  50.     RayHt,
  51.     AxisSlopeAngle,
  52.     FromIndex,
  53.     ToIndex         : BmReal;
  54.  
  55.     SpectralLine    : ARRAY [1..8] OF BmReal;
  56.     s               : ARRAY [1..MaxSurfaces, 0..5] OF BmReal;
  57.     od_sa           : ARRAY [FALSE..TRUE, 0..1] OF BmReal;
  58.  
  59.     od_fline, od_cline: BmReal;
  60.  
  61.  
  62. (** ASIN **********************************************************************
  63.  *
  64.  *  Synopsis:
  65.  *
  66.  *    PROCEDURE asin(x: BmReal): BmReal;
  67.  *
  68.  *  Description:
  69.  *    Returns the arc sine of x.
  70.  *
  71.  *  NOTE:
  72.  *    If x is out of rAnge 0.0 is returned but no further action is taken.
  73.  *)
  74.  
  75.   FUNCTION asin(x: BmReal): BmReal;
  76.   BEGIN
  77.     IF (ABS(x) > 1.0) THEN
  78.       asin := 0.0
  79.     ELSE
  80.       asin := arctan(x / sqrt(1.0-x*x));
  81.   END;
  82.  
  83.  
  84. (** COT ***********************************************************************
  85.  *
  86.  *  Synopsis:
  87.  *
  88.  *    PROCEDURE cot(x: BmReal): BmReal;
  89.  *
  90.  *  Description:
  91.  *    Resturns the cotangent of x.
  92.  *  NOTE:
  93.  *    If x is out of rAnge 0.0 is returned but no further action is taken.
  94.  *)
  95.  
  96.   FUNCTION cot(x: BmReal): BmReal;
  97.     VAR
  98.       s: BmReal;
  99.   BEGIN
  100.     s := sin(x);
  101.     IF s = 0.0 THEN
  102.       cot := 0.0
  103.     ELSE
  104.       cot := cos(x) / s;
  105.   END;
  106.  
  107.  
  108. (** TRANSITSURFACE ************************************************************
  109.  *
  110.  *  Synopsis:
  111.  *
  112.  *    PROCEDURE TransitSurface();
  113.  *
  114.  *  Description:
  115.  *    Calculate passage through a surface.  If the variable Paraxial is TRUE,
  116.  *    the trace through the surface will be done using Paraxial approximations.
  117.  *    Otherwise, the normal trigonometric trace will be done.
  118.  *
  119.  *    This routine takes the following inputs:
  120.  *
  121.  *    RadiusOfCurvature   Radius of curvature of surface being crossed.  If 0,
  122.  *                        surface is plane.
  123.  *
  124.  *    ObjDist             Distance of object focus from lens vertex.  If 0,
  125.  *                        incoming rays are parallel and the following must be
  126.  *                        specified:
  127.  *
  128.  *    RayHt               Height of ray from axis.  Only relevant if
  129.  *                        ObjDist = 0.
  130.  *
  131.  *    AxisSlopeAngle      Angle incoming ray makes with axis at intercept.
  132.  *
  133.  *    FromIndex           Refractive index of medium being left.
  134.  *
  135.  *    ToIndex             Refractive index of medium being entered.
  136.  *
  137.  *    The outputs are the following variables:
  138.  *
  139.  *    ObjDist             Distance from vertex to object focus after refraction.
  140.  *
  141.  *    AxisSlopeAngle      Angle incoming ray makes with axis at intercept
  142.  *                        after refraction.
  143.  *)
  144.  
  145.   PROCEDURE TransitSurface;
  146.     VAR
  147.       iAng,               (* Incidence angle *)
  148.       rAng,               (* Refraction angle *)
  149.       iAngSin,            (* Incidence angle sin *)
  150.       rAngSin,            (* Refraction angle sin *)
  151.       oldAxisSlopeAngle,
  152.       sagitta: BmReal;
  153.  
  154.     LABEL 999;
  155.   BEGIN
  156.     IF Paraxial THEN
  157.     BEGIN
  158.       IF RadiusOfCurvature <> 0.0 THEN
  159.       BEGIN
  160.         IF ObjDist = 0.0 THEN
  161.         BEGIN
  162.           AxisSlopeAngle := 0.0;
  163.           iAngSin        := RayHt / RadiusOfCurvature;
  164.         END ELSE
  165.           iAngSin := ((ObjDist - RadiusOfCurvature) / RadiusOfCurvature) * AxisSlopeAngle;
  166.  
  167.         rAngSin           := (FromIndex / ToIndex) * iAngSin;
  168.         oldAxisSlopeAngle := AxisSlopeAngle;
  169.         AxisSlopeAngle    := AxisSlopeAngle + iAngSin - rAngSin;
  170.  
  171.         IF ObjDist <> 0.0  THEN
  172.           RayHt := ObjDist * oldAxisSlopeAngle;
  173.  
  174.         ObjDist := RayHt / AxisSlopeAngle;
  175.         GOTO 999;
  176.       END;
  177.  
  178.       ObjDist        := ObjDist * (ToIndex / FromIndex);
  179.       AxisSlopeAngle := AxisSlopeAngle * (FromIndex / ToIndex);
  180.       GOTO 999;
  181.     END;
  182.  
  183.     IF RadiusOfCurvature <> 0.0 THEN
  184.     BEGIN
  185.       IF ObjDist = 0.0 THEN
  186.       BEGIN
  187.         AxisSlopeAngle := 0.0;
  188.         iAngSin        := RayHt / RadiusOfCurvature
  189.       END ELSE
  190.         iAngSin := ((ObjDist - RadiusOfCurvature) / RadiusOfCurvature) * sin(AxisSlopeAngle);
  191.  
  192.       iAng              := asin(iAngSin);
  193.       rAngSin           := (FromIndex / ToIndex) * iAngSin;
  194.       oldAxisSlopeAngle := AxisSlopeAngle;
  195.       AxisSlopeAngle    := AxisSlopeAngle + iAng - asin(rAngSin);
  196.       sagitta           := sin((oldAxisSlopeAngle + iAng) / 2.0);
  197.       sagitta           := 2.0 * RadiusOfCurvature * sagitta * sagitta;
  198.       ObjDist           := ((RadiusOfCurvature * sin(oldAxisSlopeAngle + iAng)) *
  199.                            cot(AxisSlopeAngle)) + sagitta;
  200.       GOTO 999;
  201.     END;
  202.  
  203.     rAng           := -asin((FromIndex / ToIndex) * sin(AxisSlopeAngle));
  204.     ObjDist        := ObjDist * ((ToIndex * cos(-rAng)) / (FromIndex *
  205.                        cos(AxisSlopeAngle)));
  206.     AxisSlopeAngle := -rAng;
  207.  
  208.   999:
  209.   END;
  210.  
  211.  
  212. (** TRACELINE *****************************************************************
  213.  *
  214.  *  Synopsis:
  215.  *
  216.  *    PROCEDURE TraceLine(line: BmInt; ray_h: BmReal);
  217.  *
  218.  *  Description:
  219.  *    Perform ray trace in specific spectral line.
  220.  *)
  221.  
  222.   PROCEDURE TraceLine(line: BmInt; ray_h: BmReal);
  223.     VAR i: BmInt;
  224.   BEGIN
  225.     ObjDist   := 0.0;
  226.     RayHt     := ray_h;
  227.     FromIndex := 1.0;
  228.  
  229.     FOR i := 1 TO CurrentSurfaces DO
  230.     BEGIN
  231.       RadiusOfCurvature := s[i][1];
  232.       ToIndex           := s[i][2];
  233.  
  234.       IF ToIndex > 1.0 THEN
  235.         ToIndex := ToIndex + ((SpectralLine[4] - SpectralLine[line]) /
  236.                    (SpectralLine[3] - SpectralLine[6])) * ((s[i][2] - 1.0) /
  237.                    s[i][3]);
  238.  
  239.       TransitSurface;
  240.       FromIndex := ToIndex;
  241.  
  242.       IF i < CurrentSurfaces THEN
  243.         ObjDist := ObjDist - s[i][4];
  244.     END
  245.   END;
  246.  
  247. BEGIN
  248.  
  249.   WriteLn('FBench Benchmark.');
  250.  
  251.   Paraxial        := FALSE;
  252.  
  253.   SpectralLine[1] := 7621.0;       (* A *)
  254.   SpectralLine[2] := 6869.955;     (* B *)
  255.   SpectralLine[3] := 6562.816;     (* C *)
  256.   SpectralLine[4] := 5895.944;     (* D *)
  257.   SpectralLine[5] := 5269.557;     (* E *)
  258.   SpectralLine[6] := 4861.344;     (* F *)
  259.   SpectralLine[7] := 4340.477;     (* G *)
  260.   SpectralLine[8] := 3968.494;     (* H *)
  261.  
  262.   (*** Load test case into working array.  The test case used in this       ***)
  263.   (*** program is the design for a 4 inch achromatic telescope objective    ***)
  264.   (*** used as the example in Wyld's classic work on ray tracing by hand,   ***)
  265.   (*** given in 'Amateur Telescope Making', Volume 3.                       ***)
  266.  
  267.   ClearAperture   := 4.0;
  268.   CurrentSurfaces := 4;
  269.  
  270.   s[1][1] :=  27.05;
  271.   s[1][2] :=  1.5137;
  272.   s[1][3] :=  63.6;
  273.   s[1][4] :=  0.52;
  274.  
  275.   s[2][1] := -16.68;
  276.   s[2][2] :=  1.0;
  277.   s[2][3] :=  0.0;
  278.   s[2][4] :=  0.138;
  279.  
  280.   s[3][1] := -16.68;
  281.   s[3][2] :=  1.6164;
  282.   s[3][3] :=  36.7;
  283.   s[3][4] :=  0.38;
  284.  
  285.   s[4][1] := -78.1;
  286.   s[4][2] :=  1.0;
  287.   s[4][3] :=  0.0;
  288.   s[4][4] :=  0.0;
  289.  
  290. (******************************************************************************)
  291. (*  Compute the looping overhead.  The Dummy procedure must have some side-   *)
  292. (*  effect so that it is not optimised out of existence.                      *)
  293. (******************************************************************************)
  294.  
  295.   StartTimer;                                   (* Start the clock.           *)
  296.  
  297.   REPEAT
  298.     Dummy;
  299.   UNTIL NullTimesUp;
  300.  
  301. (******************************************************************************)
  302. (*  Now run the benchmark.  Note that the Dummy procedure is also called so   *)
  303. (*  that we can eliminate its overhead from the looping overhead.             *)
  304. (******************************************************************************)
  305.  
  306.   StartTimer;                                   (* Start the clock.           *)
  307.  
  308.   REPEAT
  309.     (*** Do main trace in D light. ***)
  310.  
  311.     Paraxial := TRUE;
  312.     TraceLine(4, ClearAperture / 2.0);
  313.     od_sa[Paraxial][0] := ObjDist;
  314.     od_sa[Paraxial][1] := AxisSlopeAngle;
  315.  
  316.     Paraxial := FALSE;
  317.     TraceLine(4, ClearAperture / 2.0);
  318.     od_sa[Paraxial][0] := ObjDist;
  319.     od_sa[Paraxial][1] := AxisSlopeAngle;
  320.  
  321.     (*** Trace marginal ray in C. ***)
  322.  
  323.     TraceLine(3, ClearAperture / 2.0);
  324.     od_cline := ObjDist;
  325.  
  326.     (*** Trace marginal ray in F. ***)
  327.  
  328.     TraceLine(6, ClearAperture / 2.0);
  329.     od_fline   := ObjDist;
  330.  
  331.     AberLSpher := od_sa[TRUE][0] - od_sa[FALSE][0];
  332.     AberOsc    := 1.0 - (od_sa[TRUE][0] * od_sa[TRUE][1]) / (sin(od_sa[FALSE][1]) * od_sa[FALSE][0]);
  333.     AberLChrom := od_fline - od_cline;
  334.     MaxLSpher  := sin(od_sa[FALSE][1]);
  335.  
  336.     (*** D light. ***)
  337.  
  338.     MaxLSpher := 0.0000926 / (MaxLSpher * MaxLSpher);
  339.     MaxOsc    := 0.0025;
  340.     MaxLChrom := MaxLSpher;
  341.  
  342.     Dummy
  343.   UNTIL BenchTimesUp;
  344.  
  345. (******************************************************************************)
  346.  
  347.   ReportTimes;
  348.  
  349.   (*** Reference results.  These are derived from a run on Microsoft Quick  ***)
  350.   (*** BASIC on the IBM PC/AT.                                              ***)
  351.  
  352.   WriteLn;
  353.   WriteLn('Reference results:');
  354.   WriteLn('Marginal ray:           47.09479120920 0.04178472683');
  355.   WriteLn('Paraxial ray:           47.08372160249 0.04177864821');
  356.   WriteLn('Longitudinal spherical aberration:    -0.01106960671');
  357.   WriteLn('    (Maximum permissible):             0.05306749907');
  358.   WriteLn('Offense against sine condition (coma): 0.00008954761');
  359.   WriteLn('    (Maximum permissible):             0.00250000000');
  360.   WriteLn('Axial chromatic aberration:            0.00448229032');
  361.   WriteLn('    (Maximum permissible):             0.05306749907');
  362.   WriteLn;
  363.  
  364.   (*** Calculated results.  Any significant difference from the reference   ***)
  365.   (*** results should be taken very seriously.  This test has been          ***)
  366.   (*** been demonstrated to be invariant for double precision floating      ***)
  367.   (*** point formats.                                                       ***)
  368.  
  369.   WriteLn('Calculated results:');
  370.   WriteLn('Marginal ray:           ', od_sa[FALSE][0]:14, od_sa[FALSE][1]:14);
  371.   WriteLn('Paraxial ray:           ', od_sa[TRUE][0]: 14, od_sa[TRUE][1]: 14);
  372.   WriteLn('Longitudinal spherical aberration:    ', AberLSpher:14);
  373.   WriteLn('    (Maximum permissible):            ', MaxLSpher:14);
  374.   WriteLn('Offense against sine condition (coma):', AberOsc:14);
  375.   WriteLn('    (Maximum permissible):            ', MaxOsc:14);
  376.   WriteLn('Axial chromatic aberration:           ', AberLChrom:14);
  377.   WriteLn('    (Maximum permissible):            ', MaxLChrom:14);
  378. END.
  379.